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Abstract — We aim to present a tutorial on the detection, 
parameter estimation and statistical analysis of compact 
sources (far galaxies, galaxy clusters and Galactic dense 
emission regions) in cosmic microwave background ob- 
servations. The topic is of great relevance for current and 
future cosmic microwave background missions because the 
presence of compact sources in the data introduces very 
significant biases in the determination of the cosmological 
parameters that determine the energy contain, origin and 
evolution of the universe and because compact sources 
themselves provide us with important information about 
the large scale structure of the universe. 

Index Terms — 

COSMOLOGY concerns itself with the fundamental 
questions about the formation, structure and evolu- 
tion of the universe as a whole. The Cosmic Microwave 
Background (CMB) radiation is one of the foremost 
pillars of Physical Cosmology. Joint analyses of CMB 
and other astronomical observations are able to deter- 
mine with ever increasing precision the value of the 
fundamental cosmological parameters (such as the age 
of the universe, its matter and energy content and its 
geometry) and to provide us with valuable insight about 
the dynamics of the universe in evolution. 

The CMB radiation is a relic of the hot and dense 
first moments of the universe: a extraordinarily homoge- 
neous and isotropic blackbody radiation, which however 
shows small temperature anisotropics that are the key for 
understanding the conditions of the primitive universe, 
testing cosmological models and probing fundamental 
physics at the very dawn of time. For a short review on 
the CMB, see [1]. CMB observations are obtained by 
imaging of the sky at microwave wavelengths. However, 
the CMB signal is mixed with other astrophysical signals 
of both Galactic and extragalactic origin. In order to 
properly exploit the cosmological information contained 
in CMB images, they must be cleansed of these other 
astrophysical emissions first. 

Blind source separation (BSS) in the context of CMB 
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cosmology has been a very active field in the last few 
years. Most of the works have been oriented to the 
separation of the so-called diffuse components (sources 
that do not show clear borders or spatial boundaries, such 
as CMB itself and the cloud-like Galactic structures). 
A recent comparison of different BSS and other source 
separation techniques for diffuse sources can be found 
in [2]. Conversely, the term compact sources is often 
used in the CMB literature referring to spatially bounded, 
small features in the images, such as galaxies and galaxy 
clusters. For reasons that we will outline in the following 
section, compact sources and diffuse sources are usually 
treated separately in CMB image processing. We will 
devote this tutorial to the case of compact sources. As 
we will see in the following sections, and particularly in 



Section II-A many of the compact source detection tech- 
niques that are widespread in most fields of astronomy 
are not easily applicable to CMB images. In this tutorial 
paper, we will make an overview of the fundamentals 
of compact object detection theory keeping in mind at 
every moment these particularities. Throughout the paper 
we will briefly consider Bayesian object detection, model 
selection, optimal linear filtering, non-linear filtering and 
multi-frequency detection of compact sources in CMB 
images. 

I. CMB AND COMPACT SOURCES 

CMB anisotropics are extremely weak. Moreover, we 
cannot observe them directly; we observe instead a mix- 
ture of CMB and other astrophysical sources of radiation 
(usually referred to ^foregrounds or just contaminants) 
that lie along the line of sight of CMB photons. The 
foregrounds at microwave wavelengths can have Galactic 
(synchrotron radiation, brehmsstrahlung, thermal and 
electromagnetic dust emission) or extragalactic (galaxies, 
galaxy clusters) origin. Besides for ground-based exper- 
iments the atmosphere must be accounted as another 
contaminant. In addition, obviously, all observations are 
affected by instrumental noise and convolution due to 
the finite resolution of the optical devices. Therefore a 
good source separation is an indispensable prerequisite 
for any serious attempt to do CMB science. Moreover, 
the microwave window of the electromagnetic spectrum 
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has been opened for observation only very recently. As a 
result, the microwave sky is poorly known. In particular, 
the present knowledge about extragalactic sources in 
frequencies that range from 20 to ~ 1000 GHz is almost 
completely an extrapolation from observations at lower 
or higher frequencies, with only a few relevant and very 
recent direct observations [3]-[5]. 

Present-day CMB telescopes have relatively poor an- 
gular resolutions (from a few arcminutes to one degree). 
Extragalactic foregrounds show angular sizes smaller 
than, or at most comparable to, such scales^] Therefore, 
they can be considered as point-like (compact) objects 
from the point of view of CMB experiments. The de- 
tection of faint compact objects is a very important task 
that is common to many branches of astronomy. 

Compact sources constitute a special case among 
CMB foregrounds. They are the main contaminant at 
small angular scales [6] and they have a dramatic impact 
on the estimation of the cosmological parametersfrom 
the CMB signal, on the separation of other compo- 
nents and on the study of the possible non-standard 
cosmological scenarios (by probing the Gaussianity and 
isotropy of the CMB). Compact sources come in three 
main varieties: far galaxies (including radio galaxies, 
quasars, blazars, 'ordinary' galaxies, dusty galaxies with 
high stellar formation rate, protogalaxies, spheroids and 
many other kinds of objects), galaxy clusters (leaving 
an imprint on CMB radiation, mainly through inverse 
Compton scattering) and Galactic compact objects (cold 
cores, supernova remnants, dense knots of dust, etc). 
The two first classes of objects are almost uniformly 
distributed across the sky, whereas the third class is 
mainly concentrated near the Galactic plane. Practically 
all the objects belonging to the first class have angular 
sizes that are much smaller than the angular resolution 
of CMB experiments, while some galaxy clusters and 
many Galactic objects are observed as extended sources. 
The electromagnetic spectral behavior of galaxy clusters 
is pretty well known, while it is not for the other 
two classes of objects. The previous sentences serve 
to illustrate the diversity and heterogeneity of compact 
sources. Due to this heterogeneity, component separation 
techniques based on generative mixing models usually 
fail to achieve the separation of compact sources (the 
only exception to this rule are galaxy clusters). Specific 
detection (identification) and separation techniques must 
be tailored for compact sources. 

'in astronomy, the position and size of celestial objects are 
described in spherical coordinates, with center in the observer. The 
angular size or angular scale of any object is its visual diameter 
measured as an angle. Throughout this paper we will use the term 
scale as a synonym of angular scale. 




Fig. 1: A 14.6 x 14.6 square degree portion of the 
sky, centered at Galactic latitude b = 60° and Galactic 
longitude / = 100°, observed with the WMAP satellite at 
23 GHz. In the upper right corner we have superimposed 
a circle whose diameter is d = 0.88°, equal to one 
FWHM of the telescope receiver for that frequency. 

II. Basics of source detection in CMB images 
A. Particularities of the detection problem 

The detection on point-like sources is a quite old 
problem in astronomy. Over the years, many algorithms 
and software codes have been developed in order to find 
compact objects in optical images, and these methods 
have been adapted, or new others have been developed, 
every time a new section of the electromagnetic spectrum 
has become available for observation. For example, the 
widely used CLEAN algorithm [7] and its descendants 
have been at the basis of interferometric radio astronomy 
for the last 35 years. Other algorithms and codes such as 
Daofind [8] and SExtractor [9] have proven useful 
in optic, infrared and X-Ray astronomy. However, up 
to now none of these techniques (with the exception of 
SExtractor) have been much used in the context of 
compact source detection at CMB frequencies. In this 
section we will try to explain why it is so. 

Figure [T] shows a typical example of CMB image, 
taken from an observation of the sky made at 23 GHz by 
the NASA satellite WMAP [10]. At first glance it is im- 
possible to tell apart the different signals that contribute 
to the image. Besides from the detector noise, which 
operates at the pixel scale, all the other components show 
spatial correlation. In the upper right corner of Figure [T] 
we have superimposed a circle whose diameter is one 
full width at half maximum (FWHM) of the WMAP 
antenna at 23 GHz: all the features in the image (apart 
from instrument noise) smaller than those of the antenna 
are smoothed out. Moreover, that scale is the natural one 
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of point-like objects such as galaxies and galaxy clusters. 
Note that many of the fluctuations seen in the image have 
scales similar to the antenna FWHM: as we will see, this 
makes it harder to detect compact sources. 

The Clean algorithm [7] uses a simple iterative 
procedure to find the positions and strengths of these 
sources, going from the brightest intensity point in the 
image, fitting it to a source template profile, subtracting it 
from the data and iterating the whole process until some 
convergence criterion is met. In its basic implementation, 
Clean is equivalent to a least-squares minimization of 
the difference between the image and a model consist- 
ing of a sum of 'clean' source profiles [11]. Clean 
implicitly assumes that the sky can be represented by 
a small number of point sources in an otherwise empty 
field of view; unfortunately, as it is shown in FigureJT] 
this is not the case for the typical CMB image, where 
compact sources are not dominant and easily mingle with 
intrinsic CMB fluctuations of similar intensity and size. 
Even with the recent breakthroughs in the development 
of fast algorithms, CLEAN methods are computationally 
costly and require a non negligible amount of fine tuning 
work in order to perform properly on CMB images. 

Another critical aspect of compact source detection 
in CMB images is the spatial modeling of the image 
background. A common practice in astronomy is to 
subtract the image background before detecting the com- 
pact sources. It is often assumed that the background is 
smoothly varying and has a characteristic scale length 
much larger than the scale of the discrete objects being 
sought. For example, SEXTRACTOR approximates the 
background emission by a low-order polynomial, which 
is subtracted from the image [9]. Object detection is 
then performed by finding sets of connected pixels above 
some given threshold. However, again it is not the case 
for CMB images such as the one shown in FigurdT] CMB 
emission varies on a characteristic scale of order < 10 
arcminutes, comparable to the angular resolution of most 
CMB experiments, and cannot be safely modeled by a 
low order polynomial. The same applies to some Galactic 
components, especially the Galactic dust. 

The problem, in short, can be summarized as fol- 
lows: we have a unknown number of weak compact 
sources totally embedded in an omnipresent, non-sparse 
background whose intensity dominates over the compact 
sources and whose characteristic scale of variation is 
comparable to the size of the sources we are looking 
for. Under these circumstances, common methods such 
as Clean, DAOfind and SExtractor, however well 
adapted to other kind of environments may they be, find 
serious problems to work properly; we need to go back 
to the origins and review from scratch the fundamentals 



of source detection. 

As we said previously, our ability to obtain optimal 
detectors depends on our degree of knowledge about 
the probability density functions (PDFs) of signal and 
noise. Notice that in the context of compact objects in 
CMB images, the former are referred to as the signal, 
whereas the rest of the components (i.e. the CMB itself, 
the Galactic contaminants and the instrumental noise) 
are considered as a generalized noise. Indeed, only the 
brightest of the compact objects are actually considered 
as the signal, the rest of the compact object emission, that 
forms a uniform background across the celestial sphere, 
are also considered as a contributor to the generalized 
noise. Hereinafter we will refer to the generalized noise 
just as the noise, unless other case is explicitely said. 
Let us briefly describe the statistical properties of these 
two components for the case of the microwave sky: 

1) Compact sources: At the low frequencies (10— 150 
GHz) they show an almost uniform spatial distribution 
over the sky. Clustering takes an increasingly important 
role for higher frequencies. Regarding the flux distri- 
bution (number counts), the simplest models in use 
consider power law distributions, whose normalization 
and slope are determined by fitting to observational data. 
The best available models depart significantly from pure 
power law descriptions and are derived from a mixing 
of extrapolations from observations at other frequencies 
and physical modeling of the galaxies [12]. 

2) Noise: Regarding the n-point PDF, standard in- 
flationary cosmological models predict that CMB fluc- 
tuations have a multi-normal distribution^ Instrumental 
noise is also well modeled by a random Gaussian process 
(not necessarily stationary). Galactic foregrounds, how- 
ever, are strongly non-Gaussian distributed. Moreover, 
the PDF of the Galactic components is poorly known 
at microwave frequencies, particularly at small angular 
scales and, besides, it is highly non-stationary, depending 
strongly on the Galactic latitude. 



B. Two standpoints 

Three quality indicators are considered for compact 
source catalogues: reliability, completeness and accu- 
racy. Please note that astronomers do not always follow 
the standard terminology used in statistics. In fact, dif- 
ferent definitions appear in the literature. We will use the 
following terms: 

2 Some non-standard cosmological models predict some degree 
of non-Gaussianity of the CMB. However, current observational 
constrains imply that such degree of non-Gaussianity must be small. 



IEEE SIGNAL PROCESSING MAGAZINE, VOL. XX, NO. X, MONTH YEAR 



4 



1) Reliability: is one minus the ratio of false alarms 
('spurious detections') over the total number of alarms, 
for a given detection criterion: R = 1 — n e /{rid + n e ), 
where n e is the number of false alarms (spurious de- 
tections) and rid is the number of true positives. In 
some occasions, astronomers use the related term purity, 
with basically the same meaning. Reliability measures 
how many among the positives (detections) given in a 
catalogue are correctly identified as such. 

2) Completeness: is the ratio of the number of ac- 
tual detections, for a given detection criterion, over 
the true number of objects that satisfy that criterion: 
C = rid/ nt = rid/ {rid + ri m ), where n m is the number 
of true objects missed by the catalogue (false negatives) 
and therefore rit is the total number of objects that 
satisfy the selection condition. For example, we say 
that a catalogue has a completeness of 0.99 above 1 
J>0 if 99% of the sources whith fluxes above 1 Jy 
have been detected and form part of the catalogue. The 
completeness thus defined is equivalent to the sensitivity 
of a binary classification test. 

3) Accuracy: refers to the goodness of the estimation 
of the intensity -or, most commonly, flux (integrated 
intensity)- of the sources. 

The importance given to each one of these quality 
indicators depends on the research goals at each moment. 
Compact sources can be a serious nuisance for the study 
of CMB anisotropics: from the standpoint of a pure CMB 
cosmologist, the most desiderable indicator of a compact 
source catalogue will be high completeness down to the 
lowest possible flux. In a second instance, flux accuracy 
would be highly desirable. 

The opposite standpoint belongs to the astronomer 
solely interested in studying the number counts, spatial 
distribution and physical properties of objects such as 
galaxies and galaxy clusters. Possibly the astronomer 
will want to do a follow-up of the most interesting ob- 
jects at other wavelengths, then a good reliability of the 
catalogue will be the most important quality indicator, 
followed by flux accuracy and, finally, completeness. 

There is no easy compromise between the two stand- 
points. As it is well known, there is a trade-off between 
reliability and completeness. Compact source detection 
is the art of reaching the most satisfactory point of 
equilibrium between these two contradicting poles. 



3 In radio astronomy, the flux unit or Jansky (symbol Jy) is a non- 
Si unit of electromagnetic flux density equivalent to 10 -26 watts per 
square metre per hertz. 



C. Enter the detection methods: ideal optimality and 
robustness 

Up to now, we have strictly spoken about statistical 
properties of the catalogues, that is, of the output of 
the detection, without referring to the detection process 
itself. Detecting is decision making: given a image of the 
CMB, that we know for sure is corrupted by noise, how 
many (if any) compact sources are in it, and where are 
they hidden? And what are their intensities? To answer 
the first two questions, we need to choose a decision 
rule -the detector- and to somehow pre-process the 
image so that the efficiency of the detector is optimized. 
Hopefully, the way we process the image and the way we 
choose our detector are coordinated in order to obtain the 
most satisfactory results, in any of the senses described 
in the previous section. The study of this problems is 
the subject of detection theory, which is exhaustively 
covered elsewhere [13]. In this tutorial we will focus 
in the reasons why compact source separation in CMB 
observations is a particularly difficult problem, and what 
aspects of detection theory are relevant for it. 

The degree of difficulty of the detection problem is 
related to our knowledge of the signal -in this case, 
the compact sources- and the noise -everything else- 
characteristics in terms of their n-PDFs. In the ideal 
case when both PDFs are known and both of them are 
mathematically tractable it is possible, at least in theory, 
to obtain optimal detectors. When the PDFs are not com- 
pletely known, or they are not mathematically tractable, 
the determination of a good enough (not optimal any 
longer) detector is much more difficult. Assumptions 
about the signal and/or the noise PDFs are necessary to 
simplify the problem. But then, if the assumptions made 
prove to be wrong, then the outcome of the detection 
can be disastrous in terms of reliability, completeness, 
or both. On the other hand, detectors that do not make 
strong assumptions about the signal and noise PDFs are 
further away from ideal optimality, but are less sensitive 
to the accuracy of our prior knowledge about the PDFs. 
Here we find a second, empirical trade-off between ideal 
optimality and robustness. This trade-off has to do with 
the uncertainties about the signal and the noise. Note 
that it is still possible to have a detection method that 
is at the same time ideally optimal and robust, but only 
if we perfectly know the statistics of both sources and 
noise. Again, it is necessary to reach some equilibrium 
between these two opposite poles. 

III. A GUIDED TOUR FROM IDEAL OPTIMALITY TO 
ROBUSTNESS 

As we have tried to capture in the previous section, the 
current knowledge about the microwave sky is enough 
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to give us some hints about the PDF of the extragalactic 
compact sources and the noise, but there still remain 
large uncertainties about both of them. Moreover, the 
non-stationarity of the Galactic foregrounds makes it 
very difficult to devise global detection methods valid 
for all regions of the sky. Depending on how much one 
trusts the models about the signal and the noise, there 
is a panoply of different techniques that can be of use. 
In this section we will quickly overview the different 
possibilities used in the literature. 

Our starting point will be a complete description of the 
compact source detection from a Bayesian framework 
where all the information about the statistical properties 
of the sources and the background is known. Even more, 
Bayesian framework provides a unique paradigm for 
model selection. Later on, we will introduce successive 
simplifications, of the statistical model that will allow 
us to relax assumptions about the background and/or the 
signal. These simplifications will lead to loss of ideal 
optimality, but in practice to gain in robustness against 
assumptions that may be dangerously incorrect. The first 
assumption we will drop is the a priori knowledge of the 
n-PDF of the sources. This will lead to Neyman-Pearson 
detectors. Then we will relax our assumptions about the 
n-PDF of the background and try to characterize it just in 
terms of its second order statistics. This will lead to the 
well-known matched filter. Finally, we will even drop 
any pretense of knowing the covariance matrix of the 
background: we will study tools that are able to work just 
with the dispersion of the 1-PDF of the background: this 
will lead to the use of linear band-pass filters with fixed 
waveform (except for an adjustable scale parameter), for 
example wavelet functions. 

We will focus on the case of observations made at one 
single wavelength; the more complicate case of multi- 



wavelength detection will be discussed in section IV 



A. Full Bayesian approach 

Let us suppose we have a image containing a unknown 
number N of compact sources with some given spatial 
templates with compact support Tj(x), i = 1,...,JV, 
normalized for convenience to unit peak amplitude, cen- 
tered at the unknown positions Xj and having unknown 
amplitudes A\, then 



d(x) = s(x) + n(x) 



JV 

E 

i=l 



iViCx-XO + nCx), (1) 



where n(x) is a generalized noise defined as all contri- 
butions to the image aside from the discrete objects, and 
it is in general non-Gaussian distributed, non white and 
non stationary. The i th template profile can be defined in 



terms of a set of parameters 8i that include the position 
Xj of the source and some measure Ri of its size. An 
example is the circularly symmetric Gaussian profile 



t g (x; Oi) = exp 



x 



2Rf 



(2) 



This profile is a good approximation of point sources 
observed through a Gaussian beam, and in that case all 
Ri = R, the width of the beam. For the moment let us 
allow the sources to have different sizes. The whole set 
of parameters defining a single source i would be then 
@i = {Ai,9i} and the total set of unknown parameters 
belonging to the compact sources in the image would 
be = {Oi, . . . , 0jv}- Bayesian inference methods 
provide a consistent approach to the estimation of the 
parameters in a model H for the data d: 

P(d\@,H)P(@\H) 



P(@\d,H) 



(3) 



P(d\H) 

where P(Q\d,H) is the posterior probability distribu- 
tion of the parameters, P(d\Q,H) is the likelihood, 
P (Q\H) is the prior and P (d\H) = E is the Bayesian 
evidence. Inferences are usually obtained either by taking 
samples from the posterior using Markov-Chain Monte 
Carlo (MCMC) methods or by applying Maximum 
A Posteriori (MAP) inference. But for the detection 
problem we need to go further and to perform model 
selection, that is, decide between the null hypothesis 
Hq ('there are no sources present in the image') and 
an alternative hypothesis H\. This can be done by 
comparing their respective posterior probabilities given 
the observed data d: 

P{HM) = P(d|ffi)P(ffi) = E 1 P(H 1 ) 
P(H \d) P(d\H )P(H ) E P{H y K) 

where the evidences P(d\H{) = E\ and P(d\Ho) = Eq 
can laboriously be calculated as the factor required to 
normalize the posterior over 0, under the appropriate 
hypothesis. Bayesian strategy can accomplish source 
detection and parameter estimation in a single step. 
A possible strategy consists on the following iterative 
algorithm: 

1) Define the likelihood and the prior. 

2) Locate the global maximum of the posterior using 
some adequate method (MCMC sampling, down- 
hill simplex minimizer, etc). In this step it may 
be necessary to locally approximate the posterior 
distribution as a multivariate Gaussian. 

3) Estimate the parameters 0j corresponding to that 
maximum. 

4) Calculate the posterior ratio ratio ([4]) and make the 
decision between the null hypothesis (no detec- 
tion made) or the alternative hypothesis (detection 
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made). This step usually requires to approximate 
the local structure of the posterior as a multivariate 
Gaussian as well. 

5) If the detection has been made, subtract the source 
with parameters 0j from the data and do a new 
iteration from step 2. 

6) Stop the iterations when the posterior ratio crite- 
rion P(J9i|d) > P(H \d) fails. 

This procedure has some features reminiscent of the 
widely-used Clean algorithm [7], giving rise to 
Bayesian methods such as McCLEAN and MaxClean 
[14]. Obviously, the crucial points of the whole process 
are the choice of the likelihood and the prior (both for 
parameters and hypotheses), and the way the maxima of 
the posterior distribution is mapped. [14] has explored 
the problem for the case of white Gaussian noise and 
a flat prior for the parameters 0j. In a posterior work 
a much faster algorithm PowellSnakes, based on the 
simultaneous search for a number of local maxima and 
considering Gaussian correlated noise in the likelihood, 
was developed [15]. 

B. The Neyman-Pearson Lemma and Generalized Like- 
lihood Ratio Test 

The Bayesian approach is very powerful and leads 
to ideally optimal results in terms of reliability and 
completeness of the output catalogues if and only if the 
right likelihood and priors are chosen. Unfortunately, as 
we saw in section |II-A| our present knowledge of the 



statistical properties of the noise (likelihood) and the sig- 
nal (prior) is somewhat less than satisfactory. Potentially 
dangerous errors can come either from a bad modeling 
of the likelihood or from a bad modeling of the prior. 
In the current state of affairs, the uncertainties about the 
priors are larger than the uncertainties of the likelihood. 
The reason is that the assumption that the noise is 
distributed following a Gaussian with correlations is 
pretty good, at least in regions far from the Galactic 
plane where the strongly non-Gaussian foregrounds are 
subdominant. A further argument that is usually argued 
for approximating the distribution by a Gaussian grounds 
in The Central Limit Theorem, and is set as follows: 
in most of the applications, data is analysed in terms 
of their decomposition in harmonic coefficient^] rather 
than in the real space; notice that, for a given multipole 
£, the harmonic coefficient ae m is nothing but a linear 
combination of the whole observation (weighted by the 
spherical harmonic functions), and, thefore (at least for 

4 The harmonic coefficients ai m of a given signal s (x) are given 
by aim = J dfi-s (x) Yi m (a?), where the integral is defined over the 
celestial sphere, and Y( m (x) are the spherical harmonic functions. 



the most general scenarios) they are better described by 
a Gaussian distribution, than the data in the real space. 
Following this idea, in this and the following sections 
we focus on detection criteria that rely only on the noise 
statistics. The formalism below is valid for a general PDF 
of the noise, with no assumptions about Gaussianity. 

Using only information about the noise statistics, the 
Neyman-Pearson (NP) Lemma tells us that a suitable de- 
tection criterion is given by the Generalized Likelihood 
Ratio Test (GLRT): 



£(d) 



p(d|e,ffi) 
p(d|e,p Q ) 



>£*, 



(5) 



where is a constant which defines a region of ac- 
ceptance in the space of the data d and whose value is 
arbitrarily set by fixing a conveniently small probability 
of false detections due to the noise. Given the likelihoods 
P(d|0,Pi) and P(d|0,P o ) the test statistics defined 
by the GLRT Q has the maximum power for a fixed 
value of the probability of false alarms. Two possible 
ways to further increase the power of the test are to 
add further information by increasing the number of 
observable variables and to modify P(d\@,Hi) and 
P(d\@,Ho) by means of some signal processing in 
order to make it easier to discriminate between them. 
An example of this kind of procedures can be seen 
in [16], where specific linear filters were designed in 
order to optimize the power of the GLRT test for 
Gaussian noise in the space of variables (v, k), where 
v is the (normalized) intensity of the sources and k their 
(normalized) curvature, k oc — V 2 r(x) (other alternative 
definitions of the curvature of the sources can be used). 

C. Thresholding and basic linear matched filtering 

However, in many occasions it is difficult to estimate 
the distribution of quantities such as the curvature of 
the peaks of the noise (an example is the pure white 
Gaussian noise, where the curvature of the field is not 
even well defined). The most frequent situation is when 
one has to decide between the null and the alternative 
hypotheses just on the basis of pixel intensities. In that 
case, and for correlated Gaussian noise, the likelihoods 
in eq. ([5]) take the simple form 

1 



P(d|0,P o ) 
P(d|0,Pi) 



oc exp 



oc exp 



d T c M 



(d 



s) T cr 



(d-s) 



(6) 



where C is the correlation matrix of the noise, s is the 
source and the proportionality factors of both equations 
are both the same normalization, that is not essential for 
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for a fixed rate of false alarms. Figure [2] shows this effect. 
Panel (a) shows a Gaussian-shaped source embedded in 
white noise. Visually, it is almost impossible to detect 
the source. Panel (b) shows the effect of matched filter- 
ing: now the source is clearly visible among the noise 
fluctuations. Matched filtering+thresholding is generally 
considered the basic default detection technique in CMB 
astronomy [17]. The method can easily accommodate the 
possibility of having sources with different sizes/profiles, 
just by iteratively testing different r profile models and 
by looking for maxima in the signal-to-noise space [18]. 



(a) 




(b) 

Fig. 2: (a) Noisy image, (b) The same image, after a 
matched filter has been applied. 

the purpose of the GLRT test. Then, it can be proven 
that Hi has to be chosen when the statistics T(d) is 

T(d) oc d T C- l T > (7) 

for an arbitrarily chosen threshold £*. The threshold is 
usually fixed so that the probability of having a false 
alarm is conveniently small. The proportionality factor 
in this equation is a normalization that can be chosen 
arbitrarily. A usual convention is to normalize the filter 
so that the amplitude of the sources is preserved. Thresh- 
olding is the most frequently used detection technique in 
astronomy. The most typically used threshold is the so- 
called 5<T detection, that is, = 5a, where a is the 
standard deviation of the Gaussian noise. For this case, 
the statistics T(d) is the NP detector for the problem and 
the lineal operator defined by C _1 t is the well-known 
matched filter. 

The role of the matched filter in the compact source 
detection is nothing but to smooth out the noise with 
respect to the sources so that we can lower the threshold 



D. Scaling linear filters: wavelets 

Whereas the matched filter is easy to implement 
and stands clear for intuition, in many cases the noise 
correlation matrix C is not known. Under the assumption 
that sources are sparse, C can be directly estimated from 
the data d. But in some cases the sampling of d is not 
good enough to obtain a good estimation of C. A typical 
example is a image partially incomplete or covered by 
'bad' pixels. In other cases, we may suspect that C is 
not a good representation of the noise statistics but still 
want to use some kind of basic thresholding detection 
for simplicity. For those cases, wavelets (a kind of linear 
filters with a waveform that is fixed except for a scale 
parameter a) may come in handy. 

A specific wavelet basis, the Maar or Mexican Hat 
wavelet (MHW) has proved to be particularly efficient 
in many fields of astronomy, and on the detection of 
compact sources at microwave frequencies in particular. 
MHW thresholding is done after linear convolution of the 
CMB image with the MHW kernel. The thresholding is 
optimized by tuning the wavelet scale parameter a such 
that the convolved image presents a minimum standard 
deviations a, with a normalization constrain that fixes 
the observed amplitude of the source after filtering. The 
selection of the optimal value of the scale parameter a 
is driven by the statistical properties of the noise and, 
more particularly, by the typical spatial variation of the 
noise fluctuations (see for instance [19]). 

Wavelet filtering has reached a high performance 
degree within the CMB field. More efficient thresholding 
is obtained by considering an additional parameter n, that 
defines the number of oscillations of the waveform. The 
consideration of this extra parameter actually provides a 
family of wavelet basis — the MHW family [20] — , that, 
in the Fourier space, reads as: 



where n = 1 corresponds to the Fourier transform of 
the standard MHW. The additional degree of freedom 
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injected by the n parameter, allows for a better discrim- 
ination between the noise and the compact objects. 

The determination of the optimal values of the wavelet 
parameters does not depend on any previous knowledge 
of the statistical properties neither of the signal (the com- 
pact sources) nor the noise: all the relevant information 
can be directly obtained from the data d. Of course, 
different wavelets would provide different thresholding 
criteria. The optimality of the MHW family comes, on 
the one hand, from the fact that it is obtained as even 
derivatives of the Gaussian (which, we recall, reflects 
the effective compact source profile on typical CMB 
images) and, on the other hand, from empirical exercises 
that have shown how (for any practical purpose) a given 
choice of a and n provides a MHW with a shape very 
close to the matched filter [21]. 

It is worth mentioning that, as far as we are aware, 
classical wavelet analysis, like multiresolution decom- 
position, has not been used in the context of compact 
source detection in CMB images. We believe the major 
reason is related to a point previously commented: the 
generalized noise of CMB images show a lot of structure 
at both the scales of the compact sources and larger 
scales, and little is gained by applying multiresolution. 
However, multi scale approaches using wavelets have 
been applied. In particular, [19] have shown that under 
certain circumstances the characterization of the wavelet 
coefficients in a range of scales close to the optimal 
scale of detection helps to improve the photometry of 
the compact sources and the reliability of the catalogs 
of bright sources. 

Finally, it is worth mentioning that the combination 
of filtering techniques such as the matched filter or 
the MHW with astronomical detection software such 
as SEXTRACTOR [9] has very interesting possibilities. 
S Extractor gives the possibility of using a filter 
kernel defined by the user for the detection of com- 
pact sources, whereas photometry (integrated intensity 
estimation) is performed in the unfiltered map; this 
combined approach seems very promising [19], but a 
complete study of this possibility is still lacking. Note 
that in this case the improvement of the performance of 
SEXTRACTOR comes from the use of the MHW kernel. 

IV. Multi-frequency or multi-channel 

DETECTION 

Up to now, we have considered just the case of a single 
image containing a unknown number of point sources. 
Most microwave experiments, however, observe the sky 
at several different electromagnetic frequencies (chan- 
nels). In that case, the additional information coming 



from the increased number of channels may be useful to 
improve the detectability of sources. 

Two quite different cases must be considered here: 
sources whose frequency dependence is known and 
sources whose frequency dependence is unknown. The 
best example of the first case is constituted by the 
hot gas contained within galaxy clusters, that leaves 
a characteristic imprint on the CMB through inverse 
Compton scattering (thermal Sunyaev-Zel'dovich effect). 
If the frequency dependence of the compact sources is 
known, it is easy to incorporate it either in the Bayesian 
approach (see for example [14]), in the GLRT or in its 
particular case, the matched filter formalism [22]. An 
additional problem is that many galaxy clusters have 
an appreciable extension in the sky and they cannot be 
considered as strictly point-like sources. The solution is 
to introduce a free scale parameter in the detector that 
can be adjusted in order to maximize the signal-to-noise 
ratio of the detected sources. 

The problem is much more difficult when frequency 
dependence of the sources is not known. As argued 
in section II-A extragalactic sources are very hetero- 



geneous and it is impossible to define a common fre- 
quency dependence for all the many different kinds of 
objects that can be observed at microwave frequencies. 
A possibility is to give probabilistic priors for the fre- 
quency dependence in the Bayesian framework, but the 
uncertainties remain so large that this approach seems 
dangerous at this time. Very recently a non-parametric 
linear filtering scheme, based in a combination of image 
fusion and matrices of filters that are orthogonal to the 
source profiles for the off-diagonal terms of the matrix, 
has been proposed [23]. But the field of multi-frequency 
point source detection remains a largely unexplored field 
in CMB-based cosmology. 

A particular case of interest is the detection of po- 
larized point sources. The radiation coming from the 
CMB, from some of the Galactic foregrounds and from 
the extragalactic point sources is partially linearly po- 
larized; the study of CMB polarization opens a new 
window for fundamental cosmology and therefore all the 
previous discussion about the importance of detecting 
point sources serves also for the case of point source 
polarization. Given a particular coordinate system, linear 
polarization is given by the Stokes parameters Q and 
U, that are not invariant quantities but depend on the 
orientation of the detector. On the other hand, quantity 
P = \J Q 2 + U 2 is invariant and has a physical meaning. 
The detection of Q and U can be considered a particular 
case of multi-channel detection, with the particularity 
that Q and U are statistically independent quantities one 
with respect to the other. A possibility is to use any of 
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the previous techniques independently on the Q and U 
maps and then construct a 'filtered P map'. The other 
way around is to first construct a P map and then apply 
some detector to it, for example a GLRT that takes into 
account that the statistic of the noise of the P map is 
affected by the non-linear operations made to get it from 
Q and U. Both approaches are studied in detail in [24]. 

V. SOME RECENT RESULTS 

Most of the discussion we have described above has 
remained up to now almost purely theoretical, being 
the main reason that high-quality CMB measurements 
have become possible only very recently. The first all- 
sky catalogue of point sources in the frequency range 
from 23 to 94 GHz has been produced by the WMAP 
satellite [4] by means of a matched filter defined in 
harmonic space and a 5a threshold in the filtered space; 
this matched filter was defined globally for the whole 
sky. However, the statistical properties of the noise are 
highly non-uniform across the sky, due mainly to the 
presence of the Galactic foregrounds. A local approach, 
based on the MHW with index n = 2 plus a more robust 
thresholding on overlapping small areas of the sky, was 
performed in [5] and [25], leading to a new catalogue 
containing 484 radio sources whose properties in terms 
of reliability, completeness and flux accuracy have been 
determined by comparison with a complete source of 
bright objects obtained by ground-base observations at 
20 GHz (see [25] for more details). These studies have 
shown the importance of detecting locally in small areas 
of the sky instead of globally. Although the statistical and 
physical analysis of the WMAP catalogues have given 
some interesting results [26], the number of extragalatic 
sources observed so far remains small. 

VI. Some possible future developments 

CMB image processing is a relatively young area or 
research. Much work remains to be done. There is an 
incipient but resolved interest among CMB cosmologists 
to incorporate the newest ideas to solving the problem 
of compact sources in microwave images. 

Maybe the most promising idea is related to the no- 
tions of sparsity and /^-approximations. For the particular 
case of point-like objects the idea of sparse dictionaries 
comes naturally. However, the full application to CMB 
compact source detection has not been fully addressed. 
[27] have developed a methodology that minimizes the 
/p-norm with a constraint on the goodness-of-fit and 
have compared different norms against the matched filter, 
with promising results, but their simulations are not yet 
realistic enough to tell for sure the performance of the 
method under realistic circumstances. 



Another interesting possibility is to use time- 
frequency (or space-scale) representations to transform 
the data (using for example Gabor transforms or Wigner- 
Ville transforms) in order to better characterize the signa- 
ture of compact sources as compared to the generalized 
noise. Again, this line of research has just begun to be 
opened in the context of CMB image processing [28]. 

Multiscale wavelet analysis may be useful to address 
another problem that is attracting the interest of the CMB 
community: the presence of extended compact sources 
of Galactic origin in the vicinity of the Galactic Plane 
(cold cores, supernova remnants, etc). Although works 
on this topic are very preliminary, we foresee that the 
use of multiscale techniques and compact region finders 
such as S Extractor may be a good starting point to 
detect and do accurate photometry of this kind of objects. 

Finally, the field of non-linear filtering opens an al- 
most infinite range of possibilities to be explored. Some 
of the above mentioned methods, such as POWELL- 
Snakes and some of the filters used for the detection 
of polarized point sources, include non-linear filtering 
in their schemes. Non-linear filtering schemes based 
on auto-regressive models of the background PDF have 
been recently introduced in the context of diffuse source 
separation, but the application of these ideas to compact 
source detection in CMB images has not been yet fully 
addressed. We foresee that the next few years will pro- 
vide many new interesting insights, probably motivated 
by the analysis of new experiments such as the European 
satellites Planck [29] and Herschel [30] that will allow 
us to detect up to several thousands of interesting objects, 
many of them totally new. The prospective is exciting. 
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